knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix_nopushing = paste0('_sensitivityNoPushing_', as.character(iterations))
suffix_onlypushing = paste0('_sensitivityOnlyPushing_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Own Actionplan",
    'Partner Actionplan',
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,10),
  "Between-Person Effects" = c(11,17),
  "Random Effects" = c(18, 24), 
  "Additional Parameters" = c(25,25)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,10+5),
  "Between-Person Effects" = c(11+5,17+5),
  "Random Effects" = c(18+5, 24+5), 
  "Additional Parameters" = c(25+5,25+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_plan_selfPlan',
    'hu_plan_partnerPlan',
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Own actionplan",
    'Partner actionplan', 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Own actionplan",
    'Hu Partner actionplan', 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,12),
  "Conditional Between-Person Effects" = c(13,19),
  
  "Hurdle Within-Person Effects" = c(20,29),
  "Hurdle Between-Person Effects" = c(30,36),
  
  "Random Effects" = c(37, 50), 
  "Additional Parameters" = c(51,51)
  )

rows_to_pack_hu_ordinal <- list(
  "Conditional Within-Person Effects" = c(3+5,12+5),
  "Conditional Between-Person Effects" = c(13+5,19+5),
  
  "Hurdle Within-Person Effects" = c(20+5,29+5),
  "Hurdle Between-Person Effects" = c(30+5,36+5),
  
  "Random Effects" = c(37+5, 50+5), 
  "Additional Parameters" = c(51+5,51+6)
  )

WITHOUT PUSHING

Self-Reported MVPA

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_nopushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 38.46*** 2.69 [33.40, 44.04] 1.000 [0.92, 1.08] 0.000 1.003 1203 2271
Hurdle Intercept 0.25*** 0.04 [ 0.19, 0.35] 1.000 [0.84, 1.20] 0.000 1.000 1088 1791
Conditional Within-Person Effects
Daily persuasion experienced 1.02 0.03 [ 0.97, 1.08] 0.817 [0.92, 1.08] 0.973 1.002 1741 2793
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.98, 1.07] 0.837 [0.92, 1.08] 0.992 1.002 2103 3042
Daily pressure experienced 0.90* 0.04 [ 0.82, 1.00] 0.977 [0.92, 1.08] 0.302 1.000 4858 2932
Daily pressure utilized (partner’s view) 0.93 0.04 [ 0.85, 1.02] 0.940 [0.92, 1.08] 0.599 1.002 4006 2834
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.96 0.06 [ 0.85, 1.08] 0.750 [0.92, 1.08] 0.700 1.002 8120 3138
Own actionplan 1.33*** 0.06 [ 1.22, 1.46] 1.000 [0.92, 1.08] 0.000 1.000 5396 2975
Partner actionplan 1.08 0.04 [ 0.99, 1.17] 0.965 [0.92, 1.08] 0.538 1.001 6512 3309
Daily weartime NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.10 0.15 [ 0.84, 1.45] 0.758 [0.92, 1.08] 0.348 1.007 827 1540
Mean persuasion utilized (partner’s view) 1.08 0.15 [ 0.82, 1.41] 0.709 [0.92, 1.08] 0.378 1.007 833 1449
Mean pressure experienced 1.16 0.19 [ 0.82, 1.59] 0.808 [0.92, 1.08] 0.252 1.005 1105 2109
Mean pressure utilized (partner’s view) 0.96 0.16 [ 0.69, 1.34] 0.599 [0.92, 1.08] 0.346 1.005 1101 2103
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.54*** 0.11 [ 1.36, 1.79] 1.000 [0.84, 1.20] 0.000 1.000 2888 2680
Hu Daily persuasion utilized (partner’s view) 1.36*** 0.08 [ 1.21, 1.54] 1.000 [0.84, 1.20] 0.013 1.002 4621 3272
Hu Daily pressure experienced 0.97 0.13 [ 0.73, 1.28] 0.592 [0.84, 1.20] 0.792 1.000 5005 2787
Hu Daily pressure utilized (partner’s view) 1.58** 0.29 [ 1.14, 2.42] 0.997 [0.84, 1.20] 0.050 1.000 3681 1761
Hu Daily pushing experienced NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Day 0.85 0.12 [ 0.64, 1.14] 0.862 [0.84, 1.20] 0.537 1.001 7820 2875
Hu Own actionplan 8.95*** 0.86 [ 7.45, 10.88] 1.000 [0.84, 1.20] 0.000 1.000 5652 3087
Hu Partner actionplan 1.22* 0.12 [ 1.01, 1.47] 0.983 [0.84, 1.20] 0.420 1.003 5463 2654
Hu Daily weartime NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.51 0.48 [ 0.81, 2.77] 0.896 [0.84, 1.20] 0.196 1.003 1235 2156
Hu Mean persuasion utilized (partner’s view) 1.55 0.48 [ 0.84, 2.82] 0.918 [0.84, 1.20] 0.180 1.003 1187 1565
Hu Mean pressure experienced 0.35** 0.13 [ 0.17, 0.73] 0.997 [0.84, 1.20] 0.010 1.004 1541 2601
Hu Mean pressure utilized (partner’s view) 0.58 0.21 [ 0.28, 1.23] 0.923 [0.84, 1.20] 0.127 1.001 1876 2284
Hu Mean pushing experienced NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.40] NA NA NA 1.004 1351 2570
sd(Hurdle Intercept) 0.72 0.10 [0.54, 0.96] NA NA NA 1.001 1349 2021
sd(Daily persuasion experienced) 0.12 0.02 [0.08, 0.16] NA NA NA 1.001 1975 2626
sd(Daily persuasion utilized (partner’s view)) 0.08 0.02 [0.04, 0.13] NA NA NA 1.001 2475 2465
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.23] NA NA NA 1.003 1555 1609
sd(Daily pressure utilized (partner’s view)) 0.07 0.05 [0.00, 0.19] NA NA NA 1.003 1878 1845
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion experienced) 0.22 0.08 [0.07, 0.41] NA NA NA 1.003 1464 1606
sd(Hu Daily persuasion utilized (partner’s view)) 0.16 0.09 [0.01, 0.34] NA NA NA 1.003 1315 1776
sd(Hu Daily pressure experienced) 0.17 0.16 [0.01, 0.68] NA NA NA 1.000 2157 2330
sd(Hu Daily pressure utilized (partner’s view)) 0.23 0.21 [0.01, 0.89] NA NA NA 1.001 1678 1995
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.68 0.01 [0.66, 0.70] NA NA NA 1.003 6538 2471
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.78), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.71), b_hu_pressure_self_cb and
##   b_hu_persuasion_self_cb (r = 0.76), b_hu_pressure_partner_cb and
##   b_hu_persuasion_partner_cb (r = 0.75). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 110.54*** 6.14 [99.67, 123.17] 1.000 [0.94, 1.07] 0.000 1.002 744 1473
Within-Person Effects
Daily persuasion experienced 1.03 0.01 [ 1.00, 1.06] 0.972 [0.94, 1.07] 0.992 1.002 2381 2780
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.05] 0.883 [0.94, 1.07] 0.999 1.000 3188 3051
Daily pressure experienced 0.95 0.03 [ 0.89, 1.02] 0.934 [0.94, 1.07] 0.654 1.000 4197 2748
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.92, 1.05] 0.690 [0.94, 1.07] 0.927 1.001 4373 2658
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.97 0.03 [ 0.91, 1.04] 0.797 [0.94, 1.07] 0.838 1.001 5395 2627
Own Actionplan 1.07** 0.03 [ 1.02, 1.12] 0.996 [0.94, 1.07] 0.498 1.001 4381 2969
Partner Actionplan 1.05* 0.03 [ 1.01, 1.11] 0.988 [0.94, 1.07] 0.680 1.002 5206 3137
Between-Person Effects
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 1.000 3903 2616
Mean persuasion experienced 1.09 0.13 [ 0.85, 1.41] 0.763 [0.94, 1.07] 0.304 1.003 765 959
Mean persuasion utilized (partner’s view) 1.02 0.13 [ 0.79, 1.30] 0.566 [0.94, 1.07] 0.382 1.003 751 1151
Mean pressure experienced 0.91 0.13 [ 0.69, 1.22] 0.738 [0.94, 1.07] 0.285 1.002 949 1514
Mean pressure utilized (partner’s view) 1.04 0.14 [ 0.80, 1.37] 0.618 [0.94, 1.07] 0.349 1.003 873 1354
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Random Effects
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.912 [0.94, 1.07] 1.000 1.002 4283 3057
sd(Intercept) 0.29 0.04 [0.23, 0.38] NA NA NA 1.001 908 1500
sd(Daily persuasion experienced) 0.05 0.01 [0.02, 0.08] NA NA NA 1.000 1601 1315
sd(Daily persuasion utilized (partner’s view)) 0.05 0.01 [0.02, 0.09] NA NA NA 1.001 1940 2132
sd(Daily pressure experienced) 0.04 0.04 [0.00, 0.14] NA NA NA 1.004 1369 1870
sd(Daily pressure utilized (partner’s view)) 0.04 0.03 [0.00, 0.11] NA NA NA 1.000 2374 2184
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
Additional Parameters
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sigma 0.58 0.01 [0.56, 0.59] NA NA NA 1.000 5874 2881
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.88), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.84), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.78), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.87). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 3.67*** 0.11 [ 3.46, 3.89] 1.000 [-0.11, 0.11] 0.000 1.009 419 981
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.04, 0.05] 0.560 [-0.11, 0.11] 1.000 1.001 3858 3241
Daily persuasion utilized (partner’s view) 0.03 0.02 [-0.01, 0.08] 0.909 [-0.11, 0.11] 0.997 1.000 2990 2740
Daily pressure experienced -0.03 0.05 [-0.14, 0.07] 0.725 [-0.11, 0.11] 0.943 1.003 4499 2684
Daily pressure utilized (partner’s view) 0.00 0.05 [-0.11, 0.10] 0.522 [-0.11, 0.11] 0.967 1.001 4772 3218
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 0.25*** 0.06 [ 0.14, 0.36] 1.000 [-0.11, 0.11] 0.007 1.000 7438 2668
Own Actionplan 0.11** 0.04 [ 0.04, 0.18] 0.998 [-0.11, 0.11] 0.569 1.001 5645 3193
Partner Actionplan -0.01 0.04 [-0.08, 0.07] 0.584 [-0.11, 0.11] 0.997 1.000 6635 3351
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 0.40 0.27 [-0.10, 0.94] 0.943 [-0.11, 0.11] 0.099 1.007 359 598
Mean persuasion utilized (partner’s view) 0.33 0.27 [-0.19, 0.85] 0.897 [-0.11, 0.11] 0.160 1.007 362 497
Mean pressure experienced -0.34 0.29 [-0.91, 0.21] 0.887 [-0.11, 0.11] 0.158 1.007 430 640
Mean pressure utilized (partner’s view) -0.28 0.29 [-0.84, 0.28] 0.840 [-0.11, 0.11] 0.206 1.007 407 616
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.59 0.08 [0.47, 0.78] NA NA NA 1.004 602 1150
sd(Daily persuasion experienced) 0.05 0.03 [0.00, 0.11] NA NA NA 1.003 982 1807
sd(Daily persuasion utilized (partner’s view)) 0.08 0.03 [0.02, 0.13] NA NA NA 1.007 1022 681
sd(Daily pressure experienced) 0.06 0.06 [0.00, 0.24] NA NA NA 1.002 1830 2395
sd(Daily pressure utilized (partner’s view)) 0.07 0.06 [0.00, 0.24] NA NA NA 1.000 1764 1629
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
Additional Parameters
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sigma 0.96 0.01 [0.94, 0.98] NA NA NA 1.000 7952 2848
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.87), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.87), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA
Intercept[1] 3.46*** 0.99 [ 1.99, 6.09] 1.000 [0.84, 1.20] 0.000 1.001 4277 3245
Intercept[2] 7.31*** 2.21 [ 4.11, 13.20] 1.000 [0.84, 1.20] 0.000 1.001 4239 3268
Intercept[3] 19.71*** 6.39 [ 10.86, 37.64] 1.000 [0.84, 1.20] 0.000 1.000 4388 3235
Intercept[4] 85.54*** 32.07 [ 40.85, 177.13] 1.000 [0.84, 1.20] 0.000 1.001 4946 3129
Intercept[5] 2773.94*** 1761.22 [844.52, 11920.75] 1.000 [0.84, 1.20] 0.000 1.002 7124 3328
Within-Person Effects
Daily persuasion experienced 0.85* 0.07 [ 0.72, 0.99] 0.980 [0.84, 1.20] 0.580 1.001 4159 3362
Daily persuasion utilized (partner’s view) 1.01 0.09 [ 0.85, 1.21] 0.564 [0.84, 1.20] 0.948 1.000 3583 2527
Daily pressure experienced 2.02** 0.37 [ 1.34, 2.84] 0.999 [0.84, 1.20] 0.009 1.003 2038 2328
Daily pressure utilized (partner’s view) 1.13 0.24 [ 0.69, 1.77] 0.718 [0.84, 1.20] 0.527 1.002 2561 2272
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 1.64 0.53 [ 0.88, 3.08] 0.934 [0.84, 1.20] 0.148 1.000 5394 2906
Own Actionplan 1.12 0.26 [ 0.71, 1.77] 0.688 [0.84, 1.20] 0.504 1.001 4935 3695
Partner Actionplan 0.86 0.19 [ 0.55, 1.36] 0.745 [0.84, 1.20] 0.476 1.002 5147 3304
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 1.11 0.53 [ 0.44, 2.92] 0.582 [0.84, 1.20] 0.278 1.000 1543 2308
Mean persuasion utilized (partner’s view) 0.80 0.42 [ 0.29, 2.19] 0.673 [0.84, 1.20] 0.251 1.002 1841 2306
Mean pressure experienced 3.27* 1.69 [ 1.20, 9.52] 0.991 [0.84, 1.20] 0.020 1.000 1663 2685
Mean pressure utilized (partner’s view) 1.21 0.69 [ 0.41, 3.57] 0.627 [0.84, 1.20] 0.232 1.001 1989 2288
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.78 0.19 [0.46, 1.26] NA NA NA 1.000 1675 2397
sd(Daily persuasion experienced) 0.17 0.13 [0.01, 0.44] NA NA NA 1.005 785 1414
sd(Daily persuasion utilized (partner’s view)) 0.19 0.12 [0.01, 0.47] NA NA NA 1.001 1330 1750
sd(Daily pressure experienced) 0.47 0.23 [0.07, 1.01] NA NA NA 1.005 985 890
sd(Daily pressure utilized (partner’s view)) 0.30 0.29 [0.01, 1.37] NA NA NA 1.001 1106 2024
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
Additional Parameters
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.79), b_Intercept[4] and b_Intercept[3] (r = 0.86),
##   b_pressure_self_cb and b_persuasion_self_cb (r = 0.79),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.84). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix_nopushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning: There were 4 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 0.31*** 0.10 [0.17, 0.57] 1.000 [0.83, 1.20] 0.000 1.001 3635 3242
Within-Person Effects
Daily persuasion experienced 0.85 0.08 [0.70, 1.02] 0.962 [0.83, 1.20] 0.605 1.002 3818 2660
Daily persuasion utilized (partner’s view) 1.11 0.14 [0.87, 1.49] 0.797 [0.83, 1.20] 0.706 1.000 3184 2493
Daily pressure experienced 2.06* 0.53 [1.20, 4.07] 0.992 [0.83, 1.20] 0.021 1.000 2896 2535
Daily pressure utilized (partner’s view) 1.31 0.48 [0.62, 3.89] 0.775 [0.83, 1.20] 0.290 1.003 1265 560
Daily pushing experienced NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Day 1.87 0.68 [0.91, 3.78] 0.957 [0.83, 1.20] 0.094 1.003 5369 3275
Own Actionplan 1.18 0.33 [0.66, 2.05] 0.718 [0.83, 1.20] 0.418 1.001 4143 1857
Partner Actionplan 0.87 0.23 [0.52, 1.47] 0.696 [0.83, 1.20] 0.450 1.000 4581 3096
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 1.89 1.10 [0.63, 6.01] 0.870 [0.83, 1.20] 0.140 1.002 1752 2394
Mean persuasion utilized (partner’s view) 1.04 0.63 [0.33, 3.68] 0.524 [0.83, 1.20] 0.243 1.001 1993 2757
Mean pressure experienced 18.30*** 16.32 [3.32, 127.88] 1.000 [0.83, 1.20] 0.001 1.000 2763 2477
Mean pressure utilized (partner’s view) 1.13 1.09 [0.17, 7.51] 0.555 [0.83, 1.20] 0.152 1.001 2116 2956
Mean pushing experienced NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 1.21 0.25 [0.80, 1.79] NA NA NA 1.003 1504 1765
sd(Daily persuasion experienced) 0.19 0.13 [0.01, 0.49] NA NA NA 1.007 987 1080
sd(Daily persuasion utilized (partner’s view)) 0.44 0.17 [0.13, 0.87] NA NA NA 1.005 1425 1011
sd(Daily pressure experienced) 0.74 0.51 [0.04, 2.04] NA NA NA 1.002 839 1539
sd(Daily pressure utilized (partner’s view)) 0.72 0.63 [0.05, 2.65] NA NA NA 1.002 1040 969
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA
Additional Parameters
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sigma NA NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

## Warning: There were 4 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", paste0("AllModels", suffix_nopushing, ".xlsx")), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 38.46*** [33.40, 44.04] 1.000 110.54*** [99.67, 123.17] 1.000 3.67*** [ 3.46, 3.89] 1.000 NA NA NA 0.31*** [0.17, 0.57] 1.000
Hurdle Intercept 0.25*** [ 0.19, 0.35] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.02 [ 0.97, 1.08] 0.817 1.03 [ 1.00, 1.06] 0.972 0.00 [-0.04, 0.05] 0.560 0.85* [ 0.72, 0.99] 0.980 0.85 [0.70, 1.02] 0.962
Daily persuasion utilized (partner’s view) 1.02 [ 0.98, 1.07] 0.837 1.02 [ 0.99, 1.05] 0.883 0.03 [-0.01, 0.08] 0.909 1.01 [ 0.85, 1.21] 0.564 1.11 [0.87, 1.49] 0.797
Daily pressure experienced 0.90* [ 0.82, 1.00] 0.977 0.95 [ 0.89, 1.02] 0.934 -0.03 [-0.14, 0.07] 0.725 2.02** [ 1.34, 2.84] 0.999 2.06* [1.20, 4.07] 0.992
Daily pressure utilized (partner’s view) 0.93 [ 0.85, 1.02] 0.940 0.98 [ 0.92, 1.05] 0.690 0.00 [-0.11, 0.10] 0.522 1.13 [ 0.69, 1.77] 0.718 1.31 [0.62, 3.89] 0.775
Daily pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Day 0.96 [ 0.85, 1.08] 0.750 0.97 [ 0.91, 1.04] 0.797 0.25*** [ 0.14, 0.36] 1.000 1.64 [ 0.88, 3.08] 0.934 1.87 [0.91, 3.78] 0.957
Own actionplan 1.33*** [ 1.22, 1.46] 1.000 1.07** [ 1.02, 1.12] 0.996 0.11** [ 0.04, 0.18] 0.998 1.12 [ 0.71, 1.77] 0.688 1.18 [0.66, 2.05] 0.718
Partner actionplan 1.08 [ 0.99, 1.17] 0.965 1.05* [ 1.01, 1.11] 0.988 -0.01 [-0.08, 0.07] 0.584 0.86 [ 0.55, 1.36] 0.745 0.87 [0.52, 1.47] 0.696
Daily weartime NA NA NA 1.00*** [ 1.00, 1.00] 1.000 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.10 [ 0.84, 1.45] 0.758 1.09 [ 0.85, 1.41] 0.763 0.40 [-0.10, 0.94] 0.943 1.11 [ 0.44, 2.92] 0.582 1.89 [0.63, 6.01] 0.870
Mean persuasion utilized (partner’s view) 1.08 [ 0.82, 1.41] 0.709 1.02 [ 0.79, 1.30] 0.566 0.33 [-0.19, 0.85] 0.897 0.80 [ 0.29, 2.19] 0.673 1.04 [0.33, 3.68] 0.524
Mean pressure experienced 1.16 [ 0.82, 1.59] 0.808 0.91 [ 0.69, 1.22] 0.738 -0.34 [-0.91, 0.21] 0.887 3.27* [ 1.20, 9.52] 0.991 18.30*** [3.32, 127.88] 1.000
Mean pressure utilized (partner’s view) 0.96 [ 0.69, 1.34] 0.599 1.04 [ 0.80, 1.37] 0.618 -0.28 [-0.84, 0.28] 0.840 1.21 [ 0.41, 3.57] 0.627 1.13 [0.17, 7.51] 0.555
Mean pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.912 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.54*** [ 1.36, 1.79] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.36*** [ 1.21, 1.54] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 0.97 [ 0.73, 1.28] 0.592 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.58** [ 1.14, 2.42] 0.997 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.85 [ 0.64, 1.14] 0.862 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Own actionplan 8.95*** [ 7.45, 10.88] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Partner actionplan 1.22* [ 1.01, 1.47] 0.983 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.51 [ 0.81, 2.77] 0.896 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.55 [ 0.84, 2.82] 0.918 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.35** [ 0.17, 0.73] 0.997 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.58 [ 0.28, 1.23] 0.923 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 [0.23, 0.40] NA 0.29 [0.23, 0.38] NA 0.59 [0.47, 0.78] NA 0.78 [0.46, 1.26] NA 1.21 [0.80, 1.79] NA
sd(Hurdle Intercept) 0.72 [0.54, 0.96] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.12 [0.08, 0.16] NA 0.05 [0.02, 0.08] NA 0.05 [0.00, 0.11] NA 0.17 [0.01, 0.44] NA 0.19 [0.01, 0.49] NA
sd(Daily persuasion utilized (partner’s view)) 0.08 [0.04, 0.13] NA 0.05 [0.02, 0.09] NA 0.08 [0.02, 0.13] NA 0.19 [0.01, 0.47] NA 0.44 [0.13, 0.87] NA
sd(Daily pressure experienced) 0.07 [0.00, 0.23] NA 0.04 [0.00, 0.14] NA 0.06 [0.00, 0.24] NA 0.47 [0.07, 1.01] NA 0.74 [0.04, 2.04] NA
sd(Daily pressure utilized (partner’s view)) 0.07 [0.00, 0.19] NA 0.04 [0.00, 0.11] NA 0.07 [0.00, 0.24] NA 0.30 [0.01, 1.37] NA 0.72 [0.05, 2.65] NA
sd(Daily pushing experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion experienced) 0.22 [0.07, 0.41] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.16 [0.01, 0.34] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.17 [0.01, 0.68] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.23 [0.01, 0.89] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.68 [0.66, 0.70] NA 0.58 [0.56, 0.59] NA 0.96 [0.94, 0.98] NA NA NA NA NA NA NA

ONLY PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_onlypushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 57.28*** 4.53 [49.10, 66.97] 1.000 [0.93, 1.08] 0.000 1.002 1327 2040
Hurdle Intercept 3.79*** 0.98 [ 2.27, 6.26] 1.000 [0.84, 1.20] 0.000 1.005 1088 2463
Conditional Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.97 0.03 [ 0.91, 1.03] 0.856 [0.93, 1.08] 0.931 1.002 4427 2662
Daily pushing utilized (partner’s view) 0.97 0.03 [ 0.90, 1.03] 0.847 [0.93, 1.08] 0.907 1.001 3679 2873
Day 0.91 0.07 [ 0.79, 1.06] 0.894 [0.93, 1.08] 0.398 1.000 8744 2600
Own actionplan NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 1.32* 0.18 [ 1.00, 1.71] 0.977 [0.93, 1.08] 0.060 1.005 1217 2176
Mean pushing utilized (partner’s view) 1.18 0.15 [ 0.90, 1.53] 0.883 [0.93, 1.08] 0.222 1.005 1265 2076
Mean weartime NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 1.55** 0.23 [ 1.17, 2.24] 0.998 [0.84, 1.20] 0.035 1.002 3127 2531
Hu Daily pushing utilized (partner’s view) 1.84*** 0.30 [ 1.38, 2.68] 1.000 [0.84, 1.20] 0.002 1.001 2766 2989
Hu Day 0.61 0.15 [ 0.37, 1.00] 0.974 [0.84, 1.20] 0.100 1.003 7632 2768
Hu Own actionplan NA NA NA NA NA NA NA NA NA
Hu Partner actionplan NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.46** 0.14 [ 0.26, 0.81] 0.997 [0.84, 1.20] 0.018 1.001 2675 2824
Hu Mean pushing utilized (partner’s view) 0.62 0.18 [ 0.33, 1.08] 0.952 [0.84, 1.20] 0.136 1.001 2568 2699
Hu Mean weartime NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.37 0.06 [0.29, 0.50] NA NA NA 1.002 1283 2370
sd(Hurdle Intercept) 1.25 0.19 [0.93, 1.74] NA NA NA 1.002 1448 2426
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.06 0.04 [0.00, 0.16] NA NA NA 1.001 870 1486
sd(Daily pushing utilized (partner’s view)) 0.09 0.04 [0.01, 0.18] NA NA NA 1.004 977 1142
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.46 0.24 [0.07, 1.02] NA NA NA 1.000 947 1689
sd(Hu Daily pushing utilized (partner’s view)) 0.51 0.24 [0.10, 1.13] NA NA NA 1.001 1088 1582
Additional Parameters
sigma 0.67 0.02 [0.64, 0.71] NA NA NA 1.002 6374 2773
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
   
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 124.24*** 7.52 [109.69, 139.58] 1.000 [0.94, 1.06] 0.000 1.002 967 1634
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.03 0.03 [ 0.97, 1.08] 0.807 [0.94, 1.06] 0.914 1.001 3221 2733
Daily pushing utilized (partner’s view) 1.02 0.02 [ 0.97, 1.06] 0.797 [0.94, 1.06] 0.979 1.000 5065 2748
Day 0.97 0.05 [ 0.87, 1.08] 0.718 [0.94, 1.06] 0.654 1.000 7350 3336
Own Actionplan NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA
Between-Person Effects
Daily weartime 1.00** 0.00 [ 1.00, 1.00] 0.998 [0.94, 1.06] 1.000 1.001 4680 2898
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 0.97 0.06 [ 0.86, 1.10] 0.680 [0.94, 1.06] 0.635 1.001 2223 2645
Mean pushing utilized (partner’s view) 1.04 0.07 [ 0.92, 1.18] 0.739 [0.94, 1.06] 0.578 1.002 2058 2692
Random Effects
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.812 [0.94, 1.06] 1.000 1.000 3879 3246
sd(Intercept) 0.31 0.04 [0.24, 0.42] NA NA NA 1.005 890 1790
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.09 0.04 [0.02, 0.17] NA NA NA 1.004 887 821
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.03 0.03 [0.00, 0.10] NA NA NA 1.001 1398 1731
sigma 0.54 0.01 [0.52, 0.57] NA NA NA 1.001 6040 3037
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
summarize_brms(
  mood_gauss, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 3.78*** 0.12 [ 3.56, 4.00] 1.000 [-0.11, 0.11] 0.000 1.013 560 1111
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.02 0.04 [-0.05, 0.08] 0.670 [-0.11, 0.11] 0.995 1.002 2823 2591
Daily pushing utilized (partner’s view) 0.09* 0.04 [ 0.01, 0.17] 0.984 [-0.11, 0.11] 0.715 1.002 1597 2207
Day 0.28*** 0.09 [ 0.11, 0.46] 1.000 [-0.11, 0.11] 0.023 1.001 3871 3088
Own Actionplan NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 0.12 0.11 [-0.08, 0.32] 0.880 [-0.11, 0.11] 0.433 1.003 1269 1954
Mean pushing utilized (partner’s view) 0.08 0.11 [-0.12, 0.30] 0.786 [-0.11, 0.11] 0.554 1.002 1216 1958
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.61 0.08 [0.48, 0.79] NA NA NA 1.011 648 1448
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.07 0.04 [0.00, 0.16] NA NA NA 1.000 1409 1126
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.14 0.04 [0.06, 0.24] NA NA NA 1.002 1329 828
sigma 0.91 0.02 [0.87, 0.94] NA NA NA 1.001 3887 2536
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
summarize_brms(
  reactance_ordinal, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA
Intercept[1] 5.44*** 1.57 [ 3.06, 10.06] 1.000 [0.84, 1.20] 0.000 1.001 3521 2882
Intercept[2] 9.89*** 3.14 [ 5.45, 18.74] 1.000 [0.84, 1.20] 0.000 1.002 3598 2951
Intercept[3] 23.10*** 7.89 [ 12.38, 46.79] 1.000 [0.84, 1.20] 0.000 1.001 3933 2899
Intercept[4] 81.88*** 35.09 [ 38.53, 199.63] 1.000 [0.84, 1.20] 0.000 1.000 4569 3054
Intercept[5] 771.25*** 625.85 [202.52, 4855.18] 1.000 [0.84, 1.20] 0.000 1.001 5246 2545
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.31* 0.14 [ 1.04, 1.62] 0.987 [0.84, 1.20] 0.214 1.001 3641 2884
Daily pushing utilized (partner’s view) 0.96 0.13 [ 0.75, 1.25] 0.615 [0.84, 1.20] 0.810 1.002 4044 2794
Day 1.60 0.73 [ 0.68, 3.95] 0.854 [0.84, 1.20] 0.182 1.002 4822 3245
Own Actionplan NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 2.32*** 0.65 [ 1.35, 3.99] 1.000 [0.84, 1.20] 0.005 1.001 3394 2987
Mean pushing utilized (partner’s view) 1.25 0.35 [ 0.70, 2.16] 0.783 [0.84, 1.20] 0.365 1.001 3738 2562
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.85 0.21 [0.49, 1.34] NA NA NA 1.001 1569 2204
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.17 0.14 [0.01, 0.50] NA NA NA 1.004 1241 1651
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.16 0.14 [0.01, 0.57] NA NA NA 1.003 1626 1714
sigma NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
##   to inappropriate results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix_onlypushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
summarize_brms(
  is_reactance, 
  stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE Rhat Bulk_ESS Tail_ESS
Intercept 0.40*** 0.11 [0.23, 0.68] 1.000 [0.83, 1.20] 0.002 1.000 3095 2830
Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Daily pushing experienced 1.37* 0.18 [1.06, 1.87] 0.993 [0.83, 1.20] 0.154 1.000 3238 2410
Daily pushing utilized (partner’s view) 1.16 0.26 [0.80, 2.12] 0.769 [0.83, 1.20] 0.509 1.000 2310 2170
Day 1.78 0.89 [0.69, 4.81] 0.874 [0.83, 1.20] 0.160 1.000 6984 2927
Own Actionplan NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA
Mean persuasion experienced NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA
Mean pushing experienced 5.39*** 2.58 [2.38, 14.94] 1.000 [0.83, 1.20] 0.000 1.001 2293 2628
Mean pushing utilized (partner’s view) 2.37 1.15 [0.97, 6.82] 0.972 [0.83, 1.20] 0.057 1.001 2556 2640
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA
sd(Intercept) 1.52 0.31 [1.03, 2.24] NA NA NA 1.001 1431 2494
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.18 0.15 [0.01, 0.56] NA NA NA 1.003 1256 1876
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.42 0.27 [0.03, 1.15] NA NA NA 1.002 1292 1879
sigma NA NA NA NA NA NA NA NA NA
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_self_cw

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
) 

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", paste0("AllModels", suffix_onlypushing, ".xlsx")), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 57.28*** [49.10, 66.97] 1.000 124.24*** [109.69, 139.58] 1.000 3.78*** [ 3.56, 4.00] 1.000 NA NA NA 0.40*** [0.23, 0.68] 1.000
Hurdle Intercept 3.79*** [ 2.27, 6.26] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily pushing experienced 0.97 [ 0.91, 1.03] 0.856 1.03 [ 0.97, 1.08] 0.807 0.02 [-0.05, 0.08] 0.670 1.31* [ 1.04, 1.62] 0.987 1.37* [1.06, 1.87] 0.993
Daily pushing utilized (partner’s view) 0.97 [ 0.90, 1.03] 0.847 1.02 [ 0.97, 1.06] 0.797 0.09* [ 0.01, 0.17] 0.984 0.96 [ 0.75, 1.25] 0.615 1.16 [0.80, 2.12] 0.769
Day 0.91 [ 0.79, 1.06] 0.894 0.97 [ 0.87, 1.08] 0.718 0.28*** [ 0.11, 0.46] 1.000 1.60 [ 0.68, 3.95] 0.854 1.78 [0.69, 4.81] 0.874
Own actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA 1.00** [ 1.00, 1.00] 0.998 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean pushing experienced 1.32* [ 1.00, 1.71] 0.977 0.97 [ 0.86, 1.10] 0.680 0.12 [-0.08, 0.32] 0.880 2.32*** [ 1.35, 3.99] 1.000 5.39*** [2.38, 14.94] 1.000
Mean pushing utilized (partner’s view) 1.18 [ 0.90, 1.53] 0.883 1.04 [ 0.92, 1.18] 0.739 0.08 [-0.12, 0.30] 0.786 1.25 [ 0.70, 2.16] 0.783 2.37 [0.97, 6.82] 0.972
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.812 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 1.55** [ 1.17, 2.24] 0.998 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 1.84*** [ 1.38, 2.68] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.61 [ 0.37, 1.00] 0.974 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Own actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Partner actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 0.46** [ 0.26, 0.81] 0.997 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 0.62 [ 0.33, 1.08] 0.952 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.37 [0.29, 0.50] NA 0.31 [0.24, 0.42] NA 0.61 [0.48, 0.79] NA 0.85 [0.49, 1.34] NA 1.52 [1.03, 2.24] NA
sd(Hurdle Intercept) 1.25 [0.93, 1.74] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pressure experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily pushing experienced) 0.06 [0.00, 0.16] NA 0.09 [0.02, 0.17] NA 0.07 [0.00, 0.16] NA 0.17 [0.01, 0.50] NA 0.18 [0.01, 0.56] NA
sd(Daily pushing utilized (partner’s view)) 0.09 [0.01, 0.18] NA 0.03 [0.00, 0.10] NA 0.14 [0.06, 0.24] NA 0.16 [0.01, 0.57] NA 0.42 [0.03, 1.15] NA
sd(Hu Daily persuasion experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.46 [0.07, 1.02] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.51 [0.10, 1.13] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.67 [0.64, 0.71] NA 0.54 [0.52, 0.57] NA 0.91 [0.87, 0.94] NA NA NA NA NA NA NA
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()